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We investigate the eigenvalue density in ensembles of large sparse Bernoulli ran¬ 
dom matrices. We demonstrate that the fraction of linear subgraphs just below the 
percolation threshold is about 95% of all finite subgraphs, and the distribution of 
linear chains is purely exponential. We analyze in detail the spectral density of en¬ 
sembles of linear subgraphs, discuss its ultrametric nature and show that near the 
spectrum boundary, the tail of the spectral density exhibits a Lifshitz singularity 
typical for Anderson localization. We also discuss an intriguing connection of the 
spectral density to the Dedekind ? 7 -function. We conjecture that ultrametricity is 
inherit to complex systems with extremal sparse statistics and argue that a number- 
theoretic ultrametricity emerges in any rare-event statistics. 


I. INTRODUCTION 


A. Some notions about ultrametricity 


The concept of ultrametricity is related to a special class of metrics. Generally, a metric 
space is a set of elements equipped by pairwise distances. The metric d{x,y) meets three 
requirements: i) non-negativity, d{x, y) > 0 for a; 7 ^ y, and d{x, y) = 0 for a; = y, ii) 
symmetry, d{x,y) = d{y,x) and hi) the triangle inequality, d{x,z) < d{x,y) +d{y,z). 
Ultrametric spaces obey the strong triangular inequality, d{x,z) < ma,x{d{x, y), d{y, z)}, 
allowing only acute isosceles and equilateral triangles. 

The strong triangle inequality radically changes the properties of spaces in comparison 
with our intuitive expectations based of Euclidean geometry. In ultrametric spaces, the sum 
of distances does not exceed the largest summand, i.e. the Archimedean principle does not 
hold. Few counterintuitive properties are related to ultrametric balls: Every point inside a 
ball is its center; the radius of a ball is equal to its diameter; there are no balls that overlap 
only partially, if two balls have a common point, then one of them is completely embedded 
into the other. The last property implies that any ultrametric ball can be divided into a set 
of smaller balls, each of them can be divided into even smaller ones, etc. In other words, 
an ultrametric space resembles not a line of segments, but a branching tree of hierarchically 
nested balls (“basins”). The terminal points of a tree, i.e. the tree boundary, also obey the 
strong triangle inequality: for any two such points, A and B, the distance to the root of 
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the minimal subtree to which these points belong, determines the distance between A and 
B. As a result, ultrametric geometry, in general, hxes taxonomic (i.e. hierarchical) tree-like 
relationships between the elements. 

Ultrametricity is closely related to the p-adic numerical norm, which originally arose 
in the context of number theory [1]. The role of ultrametric spaces and p-adic numbers 
in number theory O |3] and in other branches of mathematics, such as p-adic analysis |1] 
and non-Archimedean analytic geometry [5], is growing. More recently, there have been 
numerous attempts to apply p-adic numbers and ultrametric spaces to various branches of 
physics ranging from string theory to inflationary cosmology (see e.g. [MI]). 

Perhaps most naturally ultrametric spaces emerge in the realm of complex systems. In¬ 
deed, the ultrametric description is inherently multi-scale, thus it is not surprising that ul¬ 
trametricity provides a fertile framework for describing complex systems characterized by a 
large number of order parameters. For instance, signatures of ultrametricity were discovered 
in the context of spin glasses (see na for an earlier review). Spin glasses are characterized 
by a huge number of metastable states and the equilibration in spin glasses proceeds via 
tree-like splitting of the phase space into hierarchically nested domains. When the temper¬ 
ature decreases, the most distant groups of spin-glass states get factorized hrst. Then, each 
of these groups splits into a number of subgroups, etc. The scales of hierarchically nested 
sets of equilibrated states satisfy the strong triangle inequality, i.e. the multi-scale splitting 
of spin-glass phases obey ultrametric relations. 

The physical origin behind the emergence of ultrametricity in spin-glass systems is quite 
general. In a complex system with an extremely large number of metastable states, the 
phase trajectories cannot totally explore all states: they cover only a tiny part of phase space 
almost without returns to already visited regions. For this reason, the factorization of spin- 
glass phases is similar to random branching in a high-dimensional (inhnitely-dimensional in 
the thermodynamic limit) space. The factorization of spin-glass phases is specihed only by 
branching points because once two branches are dispersed in a high-dimensional space, they 
never join. The ultrametric ansatz usmn proposed to describe the spin-glasses by inhnite 
(in the thermodynamic limit) number of order parameters, was proven in [T5| to be an exact 
ground state of the Sherrington-Kirkpatrick model [T6] . 

Shortly after the formulation of the ultrametric ansatz in the context of spin glasses, 
similar ideas were proposed for the description of native states of protein molecules mm- 
Proteins are highly frustrated polymer systems characterized by rugged energy landscapes of 
very high dimensionality. The ultrametricity of proteins implies a multi-scale representation 
of the protein energy landscape by means of ranging the local minima into hierarchically 
embedded basins of minima. Ultrametricity follows from the conjecture that the equilibra¬ 
tion time within the basin is signihcantly smaller than the time to escape the basin. As a 
result, the transitions between local minima obey the strong triangle inequality. In such an 
approximation, dynamics on a complex energy landscape can be modelled via a jump-like 
random process, propagating in an ultrametric space of states. The approximation of the 
protein energy landscape by a self-similarly branching tree of basins has been shown to be 
very fruitful for describing a large body of experimentally justihed features of protein con¬ 
formational dynamics and protein fluctuation mobility at temperatures covering very wide 
range from room temperatures up to the deeply frozen states [T9|. 

With respect to ultrametricity and the self-similarity of protein energy landscapes, it is 
important to keep in mind that proteins are functional systems. Proteins precisely operate 
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at the atomic level using sparse-event statistics, which is a bridge between non-local pro¬ 
tein dynamics in high-dimensional conformational space and low-dimensional movements of 
particular atomic units related to protein functionality. 


B. Ultrametricity and rare-event statistics 


In the context of high dimensionality, randomness, and rare-event statistics, ultrametric¬ 
ity emerges in diverse subjects, e.g., in genetic trees, in DNA folding (or crumpling) in 
chromosomes, in the behavior of hierarchically organized ecological, social, and economic 
systems, etc. The fact that ultrametricity is naturally rooted in high dimensionality, ran¬ 
domness, and sparse statistics, has been unambiguously observed recently in [20] • It was 
proven that in a D-dimensional Euclidean space the distances between points in a highly 
sparse sampling tends to the ultrametric distances as D —)■ cx). 

In this work we bring another intriguing connection between rare-event statistics and 
ultrametricity. We demonstrate that ultrametricity is inherit for ensembles of sparse topo¬ 
logical graphs. In particular, comparing a limiting behavior of spectral density of bi-diagonal 
random matrices, with the behavior of the function f{z) = y/--log~|r 7 (^, where ri{z) is the 
Dedekind 77 -function, and z tends to rational points on real axis, we conjecture that the 
eigenvalue distribution of sparse ensembles shares some number-theoretic properties related 
to ultrametricity. 

Our Ending provides an example of a nontrivial manifestation of rare-event Bernoullian 
noise in some collective observables in large systems (like spectral density, for example). 
From this point of view, our work should warn researchers working with efiects of ultra-low 
doses of chemicals: in the statistical analysis of, say, absorption spectra of extremely diluted 
solutions of chain-like polymer molecules, the valuable signal should be clearly purified from 
the background noise, which itself could have a very peculiar ultrametric shape. 

The paper is organized as follows. In the Section II we formulate the model under 
consideration and describe the results of our numerical simulations. In what follows we pay 
attention to statistics of linear subgraphs, presented by bi-diagonal matrices. Specifically, 
we derive the distribution of the fraction of linear subgraphs at the percolation threshold 
in the Section III, and study the spectral properties of Bernoulli ensembles of bi-diagonal 
matrices in the Section IV. The number-theoretic properties of limiting spectral statistics in 
ensembles of random bi-diagonal matrices are discussed in the Section V. The results of the 
work and the relation with other natural systems possessing the number-theoretic properties 
under consideration, are summarized in the Discussion. 


II. THE MODEL 


The Bernoulli matrix model is defined as follows. Take a large N x N symmetric matrix 
A. The matrix elements, aij = aji, are assumed to be independent identically distributed 
binary random variables such that non-diagonal elements are equal to one with probability 
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q, while diagonal elements vanish. Thus, an = 0, and for i ^ j one has 

{ 1 with probability q 

( 1 ) 

0 with probability I — q 

The matrix A can be regarded as an adjacency matrix of a random Erdos-Renyi graph, 
Q, without self-connections. The eigenvalues, A„ (n = 1, of the symmetric matrix A 

are real. Let p(A) be the eigenvalue density of the ensemble of such matrices. For iV 3> 1 
the limiting shape of p(A) is known in various cases. If g = 0(1), the function p(A) tends 
to the Wigner semicircle law, \/4:N — A^, typical for the Gaussian matrix ensembles. For 
q = c/N (c > 1), the matrix A is sparse and the density p(A) in the ensemble of sparse 
matrices has singularities at hnite values of c [2IH231- In Refs. [241127] . the behavior of the 
spectral density, p(A), has been analyzed near c = 1 (A^ 3> 1). It has been pointed out that 
the function p(A) becomes more and more singular as c approaches 1 from above. 

For ensemble of matrices A the percolation transition occurs at the value qc = 1/N. 
(The true transition occurs, certainly, in the thermodynamic limit, N —>■ oo.) For q > q^, 
the graph Q has a giant component whose size is proportional to N, and numerous small 
components (almost all of them are trees); for q < qc, the graph Q is the collection of small 
components (almost all of them are again trees). For q qc the isolated vertices dominate, 
contributing to the trivial spectral density, p(A) = 5(A). The most interesting region is close 
to qc- As long as \q — gd = the giant component is still absent. In what follows 

we discuss the behavior of the spectral density p(A) at g in this region; more precisely, we 
usually chose g slightly below qc- 

An example of a random graph with a randomly generated N x N = 500 x 500 adjacency 
matrix at g = 2.0028 x 10“^ is shown in Fig. [!} Note that the components are predominantly 
linear subgraphs. In other realizations we observed components with a single loop, but they 
are rare, usually a few per realization; more complicated components, namely those with 
more than one loop, appear very rarely. 

The whole spectrum of a symmetric matrix A in a sparse regime below the percola¬ 
tion point consists of singular peaks coming from the graph components (maximal disjoint 
subgraphs constituting the random graph). For a subgraph of size k, the correspond¬ 
ing adjacency matrix Ak of size k x k and the spectrum is determined by the equation 
det(Afc — A/fc) = 0, where Ik is the k x k identity matrix. There are k eigenvalues; apart 
from k, the eigenvalues depend on the topology of the subgraph. The heights of peaks are 
the eigenvalue multiplicities, the eigenvalue A = 0 appears with the highest multiplicity. For 
instance, the graph shown in Fig. has 222 subgraphs of size one (isolated vertices) each 
of them contributing to A = 0. There are also 10 linear chains of length three, 3 chains of 
size hve, and 1 chain of size seven contributing to A = 0; generally any chain of odd length 
has eigenvalue 0. Trees which are not linear chains also can have eigenvalue 0. One can 
bipartite any tree into two disjoint parts Xi and X 2 with edges only between Xi and X 2 . 
If the sizes ni = |Xi| and U 2 = IX 2 I of the parts are different, say ni > n 2 , the adjacency 
matrix has eigenvalue 0 with multiplicity at least ni — n 2 - 

The singular spectral density, p(A), for an ensemble of 1000 random symmetric matrices 
A, each of size 500 x 500 and generated with probability g = 2.0028 x 10“^, is shown 
in Fig. [^. The spectrum possesses the regular ultrametric hierarchical structure which 
becomes even more profound while plotting Inp(A), as shown in Fig. |^. Near the value 
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Figure 1: Typical sample of components constituting a random graph associated with an adjacency 
matrix A generated by the rule Q; the parameters are N = 500 and q = 2.0028 x 10“^. 

A = ±2 the enveloping curve of the function Inp(A) exhibits a fracture, clearly visible in 
the insert of Fig. [^, where the spectrum is replotted as p^/^(A) and the slops are depicted 
by the dashed lines Si and S 2 . (The choice of the exponent 1/3 for the curve p^/^(A) is an 
experimental selection: we have just noted that such a scaling gives the best linear shape 
for the enveloping curve and clearly demonstrates the fracture of the shape near A = ± 2 .) 

The largest eigenvalue of a tree can be roughly estimated as | = 2^/p~^^, where p is 
the maximal vertex degree of a tree (see e.g. [28H3T]). For linear chains, the maximal degree 
is p = 2 (when the chain length is at least three) and the maximal eigenvalue is = 2 . 

In the bulk of the spectrum of sparse ensemble, the randomly branching trees (with p = 3) 
and linear chains (with p = 2) dominate (some estimates are given in the following sections). 
For |A| > = 2, the linear chains do not contribute to the spectral density, and only 

graphs with p > 3 remain. For Erdds-Renyi random graphs, trees with p > 4 are rather rare. 
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Figure 2: (a) Spectral density p{X) for ensemble of 1000 symmetric random 500 x 500 Bernoulli 
matrices, generated at q = 2.0028 x 10“^ possesses regular ultrametric structure. The fracture of 
the enveloping curve near A = ±2 is clearly visible in the insert, where is plotted; (b) the 

plot of In/9(A) shows the ultrametric structure of the spectral density and terminates close to the 
maximal eigenvalue for 3-branching trees, | = 2\/2 ~ 2.83. 


so the dominant contribution comes from graphs with p = 3, giving | = 2\/2 Ri 2.83, 
consistent with simulation results (Fig. |^). 


III. DISTRIBUTION OF LINEAR CHAINS AT THE PERCOLATION 

THRESHOLD 


The spectral analysis of ensemble of random adjacency matrices A at the percolation 
threshold is cumbersome since the spectra of the components depend on their topology. For 
random Erdos-Renyi graphs, however, almost all components are trees (see e.g. [32l l33]). 
Furthermore, at the percolation point (and certainly below it) the majority of subgraphs 
(about 95% of them) are linear chains as we show below. The spectrum of each linear chain is 
easy to compute, so to get the spectral density, piin(A), of an ensemble of chain-like subgraphs 
we should know the distribution, Qn, of linear chains of length n at the percolation point. 
Below we derive an explicit expression for Qn- We hrst recall a kinetic theory approach 
which allows one to determine the size distribution of generic clusters (irrespective of their 
topology) and then adopt this framework to the computation of the size distribution of linear 
chains. 


A. General formalism: Cluster size distribution 


By dehnition, each cluster is a maximal connected component (see Fig. |^. To determine 
the cluster size distribution we employ a dynamical framework (see e.g. [33]), namely we 
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begin with a disjoint graph represented by a collection of N separated vertices and start 
linking them with a constant rate. When two vertices from distinct clusters become linked, 
the clusters join. For example, the latest link in Fig. joins two clusters of size i = 2 and 
j = 4 into a cluster of size k = i + j = 6. 


Generally, there are ix j ways to join distinct clusters of size i and j. Hence, the clusters 
undergo an aggregation process 


ihj) 


ij/(2N) 


i + j 


( 2 ) 


which proceeds with a rate proportional to the product of their sizes. We want to determine 
the cluster size distribution, Cfc(f). By dehnition, Ck{t) = Nk{t)/N, where Nkit) is the total 
number of components with k nodes. In ([^ the linking rate is chosen to be (2A^)“^. With 
this choice, the cluster size distribution Ck{t) satishes 


dck 

dt 


\ -kck 

i+j=k 


(3) 


The initial condition is Cfc(O) = 5k,i- The gain term in ([^ accounts for clusters generated 
by joining two smaller clusters whose sizes sum up to k. The second term on the right- 
hand side of Eq. (|^ represents loss due to linking of clusters of size k to other clusters. 
The corresponding gain and loss rates follow from the aggregation rule ([^. Stopping the 
aggregation process at time t leads to a random graph with adjacency matrix characterized 
by the connectivity probability q = t/N. 



Figure 3: An evolving random graph with N = 9. Links are indicated by solid lines. The graph is 
composed of three clusters: Two clusters are chains (of length 2 and 3), the third cluster is a star 
of size 4. When the link shown by a dashed line is added, the chain of length 2 and the star merge. 

Equations ([^ are well-known and they can be solved via various techniques [33] • Here we 
recall one approach as we shall use it in the next sub-section to determine the distribution 
of linear chains. The idea is to eliminate the time dependence. Solving (|^ recursively one 
gets 

ci{t) = e“*; C 2 {t) = C 3 (f) = ^ 2^3g-4t. 

This suggests to seek the solution in the form 

Ck{t) = Ckt'^-^ 


( 5 ) 
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Substituting this form into (|^, we find that the coefficients Ck can be obtained recursively 
starting from Ci = 1 and using 

(«:-l)a = t (6) 

i-\-j=k 

7 k — 2 

for k > 1. Using a generating function technique one gets Ck = Thus 

= (7) 

The large k tail is always exponential apart from the percolation point t = tc = 1 where 

Cfc(l) ~ 


B. Distribution of linear chains 


Consider now linear clusters (open chains). Let Qk{t) be the density of linear clusters of 
length k >2. We have Q 2 = C 2 , and for k > 3 the chain size distribution evolves according 
to 

2c,Qk-i + 2 Y, Q^Qj -kQk (9) 

i-\-j=k 


dt 


where the summation is taken over i > 2, j > 2. It is convenient to dehne Qi = | ci. This 
allows us to rewrite 0 as 

^ = 2Y,QiQj-kQk ( 10 ) 

i-^j=k 

here the summation runs over all i,j satisfying i + j = k. Another virtue of the agreement 
Qi = 1 Cl is that (10) is valid for all k >2. Computing the hrst few densities we get 


Qi = 5 Cl = 1 e Q 2 = C 2 = \te Qs = C 3 = e Q 3 = 1e 


( 11 ) 


etc (we have used (10) to determine Q 4 ). One is led to the conjectural behavior 


Qk{t) = lC-^e-^^ ( 12 ) 

This indeed provides the solution to and the initial condition Qa:( 0 ) = ^ 5k,i- 


Thus, the linear chain length distribution (12) is purely exponential. Rewriting (12) as 

(?>,(*) = ( 13 ) 

we see that the least-steepest decay occurs at the percolation point t = tc = 1 where 


Qk = Qk{^) “ 2 ^ ^ 


( 14 ) 
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The total chain density at the percolation point is 


+ = = 0.474928... (15) 


(15) 


Comparing this result with the total cluster density at the percolation point, ca;( 1) = 
we conclude that at the percolation point almost 95% of clusters are linear chains. This 
justihes the abundance of linear chains (see Fig. [^. 


Note that the total mass density of linear chains at the percolation point, 



(16) 


is signihcantly smaller than unity. Large clusters are usually not linear and hence despite of 
their rareness they provide a notable contribution to the mass density. 


C. Comments on more complicated components and different random graph models 

Below the percolation point, the total number of components which are not trees remains 
hnite in the thermodynamic limit —)■ cx), so it greatly fluctuates from realization to 
realization. The components which are not trees are unicyclic: each such component contains 
one loop. At the percolation point the total number of unicyclic components diverges but 
very slowly, viz. as | In A^. These results were established using probabilistic techniques, see 
[32j; they are also easily derivable using above dynamical approach (see 131 ). Overall we 
see that the contribution of the non-trees is negligible. 

The influence of the ’’shape” of the tree on the spectra is more notable. We have shown, 
however, that the majority of the clusters are not merely trees, but linear chains. In principle, 
one can modify the definition ([^ to ensure that all components are chains. If each row (and 
column) of the adjacency matrix A has no more than two non-zero elements, ^ 

for i = 1,..., A^, the corresponding graph is obviously the collection of linear chains (including 
loops) of various lengths. 

Equivalently, the chain model can be defined dynamically by using the same procedure 
as before and checking that the constraint is obeyed. In other words, we attempt to link 
vertices randomly, but we actually draw a link only if the degrees of both vertices are smaller 
than 2. The cluster size distribution characterizing the chain model evolves according to 
master equation 



(17) 


with reaction rate matrix 


/I 2 2 2 ■■■\ 



2 4 4 4 
2 4 4 4 
2 4 4 4 




•••/ 


( 18 ) 





10 


The initial condition is ca:(0) = 6 k,i- A more general model with reaction rate matrix of the 
type composed of three arbitrary rates has been studied in 


Using ([17[)-([1^ one hnds that the density of isolated vertices ci(f) and the total cluster 

(19) 


density c{t) = Cj{t) satisfy a pair of coupled rate equations: 


!| = -c.(2.-c 0. | = -i(2c-c0= 


The initial condition becomes c(0) = Ci(0) = 1. Treating c as a function of Ci we recast (19) 

c = ci\l — \ In Cl] (20) 


into Jy = ^ ~ from which 


Plugging this into the first equation (19) we obtain an implicit solution 

du 


'Cl 


v?(l — Inw) 


= t 


( 21 ) 


The following densities can be found by re-writing the rate equations (17) in a manifestly 
recursive form, 

^ = \cl- 2c2(2c- Cl), 

^ = 2CiC2 - 2c3(2c - Cl), 

= 2ciC3 -|- 2c2 — 2c4(2c — Ci), 

= 2CiC 4 -f- 4C2C3 — 2c5(2c — Cl), 

etc., and solving them recursively. 


We use the notation c^, but we emphasize that clusters generated by the above procedure 
are linear chains. There is no percolation transition in this model, so when should we stop 
the process? Recall that for Erdos-Renyi random graphs, the percolation occurs when the 
average degree is equal to one. In the model the fraction of vertices of degree zero 

(do), vertices of degree one (di), and vertices of degree two (d 2 ) are 


do — Cl 

di = 2 ^]]] Cj = 2(c — Cl) 
t>2 

d 2 = 1 — do — di = 1 — Cl — 2(c — Cl) 

The average degree (d) = di -|- 2d2 is therefore (d) = 2(1 — c). Thus, in analogy with 
Erdos-Renyi random graphs, we stop the process when the average degree reaches unity, the 
corresponding time tg is found from c{ts) = Another choice is to run the process till the 
very end when the system turns into a collection of closed loops, see [36] for this and other 
models of ring formation. 

One can generalize above model by requiring that each row (and column) of the adjacency 
matrix A has no more than p non-zero elements: X]^i — P ^ = 1,..., A^. This model 
can be investigated in the dynamic framework. Unfortunately, it has proved to be much 
more difficult for the analysis |3^ than the extreme versions corresponding to the chain 
model (p = 2) and the Erdos-Renyi variant (p = A^ — 1). 
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IV. SPECTRAL DENSITY OF CANONICAL ENSEMBLE OF BI-DIAGONAL 

RANDOM MATRICES 


A. General formalism: Spectral density 


We have seen in the Section III that the fraction of linear clusters of various lengths is 
about 95% at the percolation threshold in the ensemble of all sparse graphs. This justihes 
our study of ensembles of linear clusters (open chain and loops) only. 

It is convenient to represent the linear cluster ensemble at the percolation point (with 
the distribution Qk in chain lengths) by the bi-diagonal N x N symmetric matrix B, where 


B = 


/ 0 xi 0 0 ■ ■ ■ \ 

Xi 0 X2 0 
0 X2 0 Xs 
0 0 X3 0 


( 22 ) 


and the distribution of each xt {i = 1, is Poissonian: 


Xi = 


1 with probability q 
0 with probability 1 — q 


(23) 


Due to ergodicity, the spectral density, piin(A), of the ensemble of iV x V matrices B in 
the limit iV —)■ cx) coincides with the spectral density of an individual matrix B and can be 
computed straightforwardly. Note that for any Xk = 0, the matrix B splits into independent 
fragments. So, it is instructive to consider subchains consisting of sequential set of Xk = ■ 
The probability to have a subchain with n consecutive 1 is Qn = q'^, as it follows from (23). 
Comparing with the distribution of linear subchains in ultra-sparse matrix 

ensemble at the percolation point (see 0) , we conclude that q = e 

In what follows we derive piin(A) for arbitrary values of q. This allows us to investigate 
the spectral density of the ensemble of bi-diagonal random operators in the limit g —)• 1 and 
uncover it number-theoretic structure. 

The spectrum of linear chain of length n (i.e. of symmetric n x n bi-diagonal matrix B 
with all Xk = f, k = 1, ...,n) reads 


Afr = 2 cos 


irk 

n + 1 ’ 


(k = l,...,ij) 


(24) 


The spectral density of the ensemble of N x N random matrices B with the Poissonian 
distribution of matrix elements can be written as 


N 




, k=l 


= lim —r-Im {Gai{X — is)) 

N—¥oo 'JYJSJ 
€— >-0 

N 

= lim — ImG„(A-ie) 

N^oo 'J^ f ^ 

£->•0 n=l 


( 25 ) 
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where we have used an identity 


lim Im --^ 

X i IS 


=Fz7r5(A) 


The Green function 


Gr,{\-ze) = J2 


k=l 


A — Ai 


le 


(26) 


with Afc given by (24), is associated with each particular gapless matrix B oi n sequential 
ones on the sub-diagonals, while (...) means averaging over the distribution = q'^- 


Collecting (24), (25) and (26), we hnd an explicit expression for the density of states, 

Piin(A): 

N n ^ 




iv—^oo 'ji _/y 
£->•0 n=l 


^ (A-2cos;^) +£= 


(27) 


The sample plots piin(A) for two different values of q, namely for q = 0.7 and g = e ^ ~ 
0.37 computed numerically via (27) with £ = 3 x 10“^ are shown in the Fig. 




(b) 


Figure 4: The spectral density piin(A) for the ensemble of bi-diagonal matrices at g = 0.7 (a) and 
q = e~^ ~ 0.37 (b); the regularization parameter e is taken £ = 2 x 10“^. 


B. Analysis of enveloping curves and tails of spectral density /3iin(A) 


One can compute the enveloping curves for any monotonic sequence of peaks. In Fig. 
we show two series of sequential peaks: Si = {1-2-3-4-5-...} and S 2 = {2-6-7-...}. Any 
monotonic sequence of peaks corresponds to the set of eigenvalues Xk constructed on the 
basis of Farey sequence [38]. For example, as shown below, the peaks in the series Si are 
located at A^ = —2 cos {k = 1, 2,...), while the peaks in the series S 2 are located at A*,/ = 
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— 2 cos = 2,3,...). Positions of peaks obey the following rule: let {Xk-i, Xk, Afc+i} 

be three consecutive monotonically ordered peaks (e.g., peaks 2-3-4 in Fig. |^, and let 


Afc-i 


—2 cos 


npk-i 

Qk-i 


Afc+i — 


—2 cos 


Qk+l 


where pk and qk (fc = 1,..., N) are coprimes. The position of the intermediate peak, Xk, is 


\ o "^Pk Pk Pk—1 Pk+1 Pk—1 T Pk+1 

Xk = —2 COS - , — = - © - = - 

qk qk qk -1 qk+i qk-i + qk+i 


( 28 ) 


The sequences of coprime fractions constructed via the © addition are known as Farey 
sequences. A simple geometric model behind the Farey sequence, known as Ford circles 
[39] . is shown in Fig. |^. Take the segment [0,1] and draw two circles Oi and O 2 both of 
radius r = ^, which touch each other, and the segment at the points 0 and 1. Now inscribe 
a new circle touching Oi, O 2 and [0,1]. Where is the position of the new circle along 
the segment? The generic recursive algorithm is shown in Fig. and constitutes the Farey 
sequence construction. 



(a) 


(b) 


Figure 5: Ford circles as illustration of the Farey sequence construction: a) Each circle touches two 


neighbors (right and left) and the segment; b) Determination of the position of newly generated 
circle via the © addition: ^ ^ . 

qk-1 Qk+l Qk-l+Qk+l 


Consider two examples related to the set of eigenvalues. 

• The monotonic series Si is can be constructed recursively by taking 

TT TT 

Ai = —2 cos = 0 ) Aoo = —2 cos — = 2 


1 1 

"'2®I 


The next-to-highest peak 2 is located at 

A 2 = —2 cos 

Peak 3 is located at 

X " / 2 1 

A 3 = -2 cos ^ ( 2 © Y 


27r 

= —2 cos — = 1 
3 


37r r- 

= —2 cos — = a /2 
4 
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Peak 4 is located at 

A 4 = —2 cos 


3 1 


4:71 ^/5+l 


etc. The generic expression for Afc in the series Si reads 

kiT 


Xk = —2 cos 


k + V 


(*:= 1 . 2 ,...) 


(29) 


The monotonic series S 2 is can be constrncted similarly by taking as reference points 
Ai and A 2 . For peak 6 we have 


Afi = —2 cos 


1 2 
"'2®3 


. Stt Vb-l 

= —2 cos — = - 

5 2 


Peak 7 is located at 


A? = —2 cos 


1 3 

x,-©- 


dvr 

= —2 cos — 
7 


etc. The generic expression for Afc' in the series S 2 reads 

k'n 


Xu = —2 cos 


2k' 


{k' = 2,3,...) 


(30) 


As follows from (27), the fnnction piin(A) is nonzero only at Xk- In Fig. |^we show the set 
of eigenvalnes A„(/c) = —2 cos for n = 1 ,..., 50. Every point designates some eigenvalne, 
Xk (see (24)); the points along each curve correspond to different values k for a specific value 
n. The collection of eigenvalues in each horizontal line, coming from different n, gives the 
multiplicity of corresponding eigenvalue Xk and contributes to the height of a peak in the 
spectral density piin(A). 

The heights of peaks in any monotonic sequence (like Si or S 2 ) can be computed directly 
from (27) by summing geometric series along each horizontal line in Fig. For example, 
for the sequences Si and S 2 we have: 


Ai = —2 cos I 
Ao = —2 cos ^ 


5i: 


p'^i(Ai) = + ... 

P'^n-^2) = + g® + g^^ + ... 


(31) 


A, = -2cos^ p^fiXk) = E = lii {k = 1,2,...) 


S=1 


and 


& 


A 2 = —2 cos ^ 
As = —2 cos ^ 


Xk' = —2 cos 


nk' 

2k'-l 


= g^ + g^ + g® + g^^ + ... 
p^'^iXs) = q'^ + q^ + q^^ + g^® + ... 


s=l 


( 32 ) 


(k = 2,3,...) 
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2 = -2cos{7t) 

(V5+1)/2 = -2cos(47i/5) 

V2 = -2cos(3ii/4) 

1 = -2cos(2W3) 

(V5-1)/2 = -2cos(37i/5) 

0 = 2cos(7i/2) 

-{V5-1)/2 = 2cos(37t/5) 

-1 = 2cos(2;i/3) 

-V2 = 2cos(3;t/4) 

(-V5+1)/2 = 2cos(4;t/5) 

-2 = 2cos(jr) 

Figure 6: Family of 50 curves Wn{k) = —2cos^^. Each curve corresponds to a particular 
n = the points mark values k = Each horizontal line is the multiplicity of 

the corresponding root and contributes to the height of the peak in the spectral density. 



Note that the expressions (31 )-([3^ are tightly linked to so-called visibility diagram (our 
notion is slightly modihed with respect to a visibility diagram dehned in mi). Consider the 
square lattice of integer points (m, n) and add a weight to each vertical row. Emit rays 
at rational angles ap^q = arctan ^ from the point (0,0) within the wedge [vr/d, 7r/2] in the 
positive direction as shown in Fig. 

Let us sum up the weights, g"*, of integer points along each ray emitted at the rational 
slope 2. For example, consider the sequence of slopes tanai ^ = | with k = 1,2,... and 

tanafc/ 2 fc'-i = 2 k''-i shown in Fig. 7 One can verify that these geometric series exactly 
reproduce the sequences Si and S 2 m Eqs. (31)~(32) at corresponding k and k'. 

Two parametrically dehned enveloping curves, namely p^^{k)') for the sequence S 
{k = 1,2,...) and (A^/,p'^2(fc')) for the sequence S 2 {k' = 2,3,...), are shown in Fig. 
dashed lines. Equations (31)-(32) allow to hnd the singular tails of distributions 
A —)■ 2 and p^^{X) as A —)■ 0. Expand now the parametrically dehned curves in (31)-(32) at 
fc —)■ 00 and k' —)■ 00 : 


4ji by 
20 as 




^ = -2cos^ 


k^OQ 




So '■ X' = —2 cos 


irk' 


2k'-l 


k'^00 


— 


^^(A 


k^oo 




p 


k'; 


k'^oo 


q 


2k' 


Expressing k and k' in (33) in terms of A and A', we get 

k 


TT 


v/2^’ 


k' —)■ — 
^ ^2A' 


(33) 


(34) 


Asymptotics (34) together with (33) permit us to hnd the singular tails of the eigenvalue 
distribution piin(A) at A —)■ 2 and A' —)• 0: 

p®;(A ^ 2) ^ ^ 0) ^ 


(35) 
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1/1 ■■■ 5/6 4/5 3/4 2/3 3/5 4/7 ■■■ 1/2 



Figure 7: Visibility diagram. Each point in the vertical row m carries a weight q"^. Integer points 
within the wedge [7r/4,7r/2] are designated by circles. The dashed rays are drawn at rational angles 
Oip,q = arctan^. The weights corresponding to marked integer points, are summed up along 
the rays. 


The edge singularity —)• 2) at A —)• 2 in Eq. (35) of the eigenvalue distribution 

reproduces the corresponding Lifshitz tail of the one-dimensional Anderson localization in 


_pj-d/2 


where E = 2 — X and d = 1. The 


random Schrodinger operator IT2]: p{E) 
appearance of the edge singularity p{E) ~ g-i/V® situation is purely geometric, it is 

not relied on any entropy-energy-balance consideration like optimal fluctuation 


V. SPECTRAL DENSITY OF BI-DIAGONAL MATRICES IN THE LIMIT q 1 

AND THE DEDEKIND r?-FUNCTION 


A. Limiting properties of the Dedekind //-function via duality relation 


Recall the dehnition of the Dedekind //-function |13]: 

OO 

= z = x + iy {y > 0) 

n=0 

Consider now the following function 

f{z) = A-^\p{z)\{lmzY^^; A = 

This function satisfies the following duality relation 





1 /^ 


^ ^ 0.77230184... 


(36) 


(37) 


ms — kr = 1; {m, s, k,r} e Z; y>0 (38) 
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where (x) , {f } denote fractional parts of corresponding qnotients. Two particnlar identi¬ 
ties follow from 


/ 




/ 


k + 1 
k' 

2k' - 1 


+ w] = f 


k + 1 {k + lYy 

2 i 

+ 


2k' - 1 {2k' - Ify 


(39) 


Using (39) we can obtain the asymptotics of rj{z) when y —)• O’*". Take into acconnt the 
relation of Dedekind r^-fnnction with Jacobi elliptic fnnctions: 


^;(0,e"^) =rf{z) 


where 




du 




u=0 


X;(-ir(2n+i), 


,7rm(n+l)2 


(40) 

(41) 


Rewrite (39) as follows: 


V 


k 


k + 1 


+ iy 


VWTiyy 


n=0 


V 


+ 


fc -|- 1 {k + l^y 


(42) 


(Since the compntations are identical for both eqnations in ( pl9| ), we do not reprodnce the 
ion for tl 
, we get: 


derivation for the second one, and will write for it the hnal answer only). Applying (40)-(41) 
to 


V 


k 


k + 1 


+ iy 


1 




g4®(fe+i + (fe+i) 2 yj ^^(—i^"-(2n l)e 




n=0 


1/3 


(43) 


Eqnation (43) enables us to extract the leading asymptotics of |?7 -1- iy) | in the i/ —?• O’*" 

limit. Noting that every term of the series in (43) for any n > 1 and for small y (i.e. for 
y ^ {k + 1)~^) converges exponentially fast, we have 


k + 1 


+ iy 


-)■ 


1/-40+ Vi^ + ^)y 


g 12(k + l)'^y j 


(44) 


Similarly, we get the asymptotics for \r] + w) \ at i/ —)■ O’*" (valid for y {2k' — 1) ^). 

Thus, we have hnally: 


In 


A: -I- 1 


+ iy 


y^Q+ 


In 


k' 


2k' - 1 


+ iy 


-)■ 


-)■ 


k + lV I2y 


TT 


y-S>0+ 


2fc' - 1 Y I2y 


TT 


(45) 


In (45) we have dropped out the subleading terms at i/ —)■ 0+, which are of order of 0(lni/). 
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B. 


Tails of bi-diagonal matrices at —?> 1 and relation with Dedekind r/-fnnction 


Expressions (31)-(32) allow to analyze the tails of the distribntion piin(A) when g —)■ 1 


and k is hxed. Namely, taking the limit g —)■ 1 in (31)-(32), we arrive at: 




g" 


1 — g*^+^ 

q2k'-2 


-)■ 


q^l 


q 


2k'-l 


(fc + l)(l-g) 
1 


(46) 


-)■ 


q^l 


(2 /c'-l)(l-g) 


These expressions lead us to the following asymptotics of tails of the distribution p(A) at 

72^ 


ARi2 — —^2 and A' ^ 0: 

k^ lfc>l 2k' lfc>l 


pfM ,(A^ 2 )^ , 

^ vr(l-g) 

P? I 1 (A' 0 ) -)■ — 2 ^-r 

’ 7r(l-g) 


(47) 


Note the essential difference with (35). It is easy to see that all tails of intermediate mono¬ 
tonic sequences have scaling behavior as in the sequence pf^(A). 

Return now to (4^ and substitute for k and k' the asymptotic expansions k{X) ~ 
and k'{\') k, | (see (|34p). Since near the spectrum edges, /c S> 1 and fc S> 1, we get 


-\n\ri {\ + iy)\ 


A ^2 

y ->-0 


\/-ln|? 7 (A + M 


A 


^^0 ^ V127r2/ 


(48) 


The tails in (48) exactly match the tails of the spectral density p"^(A) and p^(A) in (47) if 


we make an identification 1 — g = y/VlTry. 

The same analysis can be carried out for any monotonic sequence of peaks of the spectral 
density depicted in Fig. This brings us to the central conjecture od this work, which 
uncovers the number theoretic nature of the spectral density of Poissonian ensemble of large 
random bi-diagonal symmetric matrices. 

Conjecture. Let piin(A, g) be the spectral density of ensemble of iV x bi-diagonal symmet¬ 
ric random Bernoulli matrices with the Poissonian distribution of nonzero matrix elements 


defined in Eqs.(22)-(23). For N ^ 00 we have the relation: 

Plin(A, g) 


lim 


-In p fAarccos(-A/2) 


= 1 


(49) 


where rj is the Dedekind //-function defined in (36). One can say, that Eq.(49) establishes the 


link between spectral density of random one-dimensional Schrodinger-type operators with 
the limiting values of modular function at rational points on the real axis. 

Empirical evidence in favor of the above conjecture is presented in Fig. where we plot 
three functions: 
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Fig. 


Numerically generated spectral density of ensemble of iV x iV {N = 3000) bi-diagonal 
Bernoulli random matrices B (see (22)) with the Poissonian distribution of “1” at 
q = 0.98; 


Fig. 1^: Numerically computed eigenvalue density p{\,q) given by Eq.(27) for matrices of size 
N = 3000 and q = 0.98; 


Fig. 1^: Analytic expression given by the right-hand side of rt49h ai y = 




12n 


10 


-5 



(a) (b) (c) 


Figure 8: Spectral density piin(A) of ensemble of bi-diagonal Poissonian random matrices (shown 
only for A > 0): a) direct numeric simulation; b) numeric visualization of Eq. (27); c) analytic 
construction provided by Eq. (49). 


VI. DISCUSSION 


A. Ultrametricity and rare-event statistics 


We investigated the spectral density of a sub-ensemble of sparse symmetric Bernoulli 
matrices near the percolation point. We showed that in the vicinity of the percolation 
point linear chains constitute the overwhelming majority of clusters—about 95%. This led 
us to study the spectral density of ensemble of linear chains only. The advantage of this 
consideration is the possibility of exact treatment which allowed us to uncover the number- 
theoretic structure of the corresponding density of states. 

In a more practical setting, our analysis is applicable to investigation of polydisperse 
solutions of linear polymers at any concentrations with Poissonian distribution in chain 
lengths, where the ultrametricity could be directly seen in the absorption spectra. One 
tacit message coming from our work is that experimenting with biological activity of highly 
diluted solutions of biologically active substances, one should pay attention to a very peculiar 
structure of background noise originated from the rare-event statistics of dissolved clusters. 
In order to make conclusion about any biological activity of regarded chemical substance, 
the signal from background noise should be clearly identihed. 
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As concerns exact results, we showed that the spectral density, piin(A), in the ensemble of 
bi-diagonal random symmetric matrices with the Poissonian distribution of matrix elements 
demonstrates a very peculiar ultrametric structure. 

The derivation of the shape of the enveloping curve for the spectral density, piin(A), of 
random Shrodinger-type operators presented by random symmetric bi-diagonal matrices, is 
based on simultaneous localization of eigenvalues, \k = 2cos| (fc = 1,2,...), constituting 

/c 

the principal sequence, and of their multiplicaticity (degeneracy), p{Xk) = (0 < g < 1). 

Eliminating k from the parametric dependence (A*,, p(Afc)), one gets the signature of a Lifshitz 
tail for Anderson localization at the edge of the spectrum, p(A —)■ 2) ~ qr-W 2 -A^ typical 
for the one-dimensional Schrodinger random operators. We emphasize that in our case the 
derivation is purely geometric. 

Analyzing the limit g —)■ 1, we have demonstrated that the function piin(A) shares number- 
theoretic properties connected with the structure of the Dedekind p—function near the real 
axis at rational points. Intriguingly, the limiting behavior of the Dedekind function near 
real axis at rational points has tantalizing connections with Vassiliev knot invariants m 
and with quantum invariants of 3-manifolds [15]; a new concept of quantum modular forms 
was recently suggested |16] as a potential explanation. 


B. Manifestation of the Dedekind 77 -fnnction in natural science 


Finally, we mention two other physical systems exhibiting a connection with the rj func¬ 
tion. In Ref. m it has been argued that buckling of a leaf of some plants (like lettuce or 
spinach) is related with the isometric embedding of hyperbolic graphs (Cayley trees) into 
the 3D Euclidean space. This buckling is described by a function, g{z)^ (Jacobian of trans¬ 
formation), which dehnes the deformation of the hyperbolic graph in its projection from the 
3D Euclidean space into the complex plane (see also the Section 5.1 of |1H])- 


For a regular 3-branching Cayley tree, the function g{z) can be expressed via the Dedekind 
g-function gziiin]. It has been noted in ITHI t hat properly normalized function g{z), namely, 
f{z) = A“^|? 7 (^)|(Im^)^/^, dehned in Eq.(37), can be regarded as a continuous tree isomet- 
rically embedded into the half-plane Im^r > 0. The constant A in (37) is chosen to set the 
maximal value of the function f{z) to 1: 0 < f{z) < 1 for any 2 ; in the upper half-plane 
Im ;2 > 0. The function f{z), being the modular function, has the following property |19]: 
f{z) has local maxima equal to 1 at the vertices of the 3-branching Cayley tree (and only of 
them) isometrically embedded into the upper half-plane Im^; > 0. The cut of the function 
f{z) G [0.985,1.0] is shown in Fig. [^, while in Fig. the same function is replotted in the 
unit disc, obtained via the conformal mapping of the upper half-plane Im 2 ; > 0 to the unit 
disc. Note the striking similarity of Fig. with the set of touching Ford circles (Fig. |^. In 
fact, the a:-coordinates of centers of all lacunas in Fig. obey the Farey construction. 


Another well-known manifestation of number theory in natural science arises in the con¬ 
text of phyllotaxis PI- The connection of the cell packing with the Farey sequences has 
been observed long time ago. It was not clear, however, why nature selects the Fibonacci 
sequence among other possible Farey sequences. A tantalizing answer to this question has 
been given by L. Levitov, see |5T]. He proposed an “energetic” approach to the phyllotaxis. 
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V5-1 

2 



(b) 


Figure 9: Cut {f{z) G [0.985,1.0]) of the contour plot of the function f{z) (see ©): a) the 
representation in the upper half-plane, z = x + iy; h) the representation in the unit disc obtained 
by the conformal mapping of the upper half-plant to the unit disc. 


suggesting that the development of a plant is connected with an effective motion along the 
geodesics on the surface associated with the energetic relief of growing plant. 

Intriguingly, the energy relief discussed in [51] coincides with the relief of the function 
f{z) depicted in Fig. The black arcs drawn in Fig. |^,b are parts of semicircles orthogonal 
to the boundaries and represent the geodesics of the surface f{z), which are open level lines 
(by the symmetry, there are only two such level lines in Fig. and three in Fig. |^). The 
best rational approximation of the boundary point, at which the solid geodesics in 

Fig. 1^ terminates, is given by the continued fraction expansion 


Vs - 1 

2 


1 


1 + 



(50) 


cut at some level k. Cutting (50) at sequential k, we get the set of fractions constituting 
the interlacing Fibonacci sequences, {FI}, in nominators and denominators: 


1 1 2 3 5 Fk-i Fk 
V 2’ 3’ 5’ 8’-’ Fk ’Ffc+i 


Thus, the normalized //-function, f{z), plays a role of the energy relief of a growing plant. 

One can say that the Fibonacci numbers are the benchmarks, which £x the connection 
between ultrametric and hyperbolic geometry, since, on one hand they appear as discrete 
symmetries of the hyperbolic space, and on the other hand, being a subset of p-adic numbers, 
they uniquely parameterize the distances in the ultrametric space [52] . 
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